Exercise 4.1¶

import gdal # Import GDAL library bindings
import numpy as np

# The file that we shall be using
# Needs to be on current directory
filename = 'files/data/MCD15A2.A2011185.h09v05.005.2011213154534.hdf'

g = gdal.Open(filename)
# g should now be a GDAL dataset, but if the file isn't found
# g will be none. Let's test this:
if g is None:
    print "Problem opening file %s!" % filename
else:
    print "File %s opened fine" % filename
    
    
subdatasets = g.GetSubDatasets()
for fname, name in subdatasets:
    print name
    print "\t", fname

# Let's create a list with the selected layer names
selected_layers = [  "Lai_1km", "FparLai_QC" ]
# We will store the data in a dictionary
# Initialise an empty dictionary
data = {}
# for convenience, we will use string substitution to create a 
# template for GDAL filenames, which we'll substitute on the fly:
file_template = 'HDF4_EOS:EOS_GRID:"%s":MOD_Grid_MOD15A2:%s'
for i, layer in enumerate ( selected_layers ):
    this_file = file_template % ( filename, layer )
    print "Opening Layer %d: %s" % (i+1, this_file )
    g = gdal.Open ( this_file )
    
    if g is None:
        raise IOError
    data[layer] = g.ReadAsArray() 
    print "\t>>> Read %s!" % layer

# scale the LAI    
lai = data['Lai_1km'] * 0.1

# pull out the QC 
qc = data['FparLai_QC']
# find bit 0
qc = qc & 1

# generate the masked array
laim = np.ma.array ( lai, mask=qc )

def read_modis_lai(filename):
    '''  Read MODIS LAI data from filename
         and return masked array
    '''

    g = gdal.Open(filename)
    # g should now be a GDAL dataset, but if the file isn't found
    # g will be none. Let's test this:
    if g is None:
        print "Problem opening file %s!" % filename
    else:
        print "File %s opened fine" % filename
        
        
    subdatasets = g.GetSubDatasets()
    #for fname, name in subdatasets:
        #print name
        #print "\t", fname
    
    # Let's create a list with the selected layer names
    selected_layers = [  "Lai_1km", "FparLai_QC" ]
    # We will store the data in a dictionary
    # Initialise an empty dictionary
    data = {}
    # for convenience, we will use string substitution to create a 
    # template for GDAL filenames, which we'll substitute on the fly:
    file_template = 'HDF4_EOS:EOS_GRID:"%s":MOD_Grid_MOD15A2:%s'
    for i, layer in enumerate ( selected_layers ):
        this_file = file_template % ( filename, layer )
        #print "Opening Layer %d: %s" % (i+1, this_file )
        g = gdal.Open ( this_file )
        
        if g is None:
            raise IOError
        data[layer] = g.ReadAsArray() 
        #print "\t>>> Read %s!" % layer
    
    # scale the LAI    
    lai = data['Lai_1km'] * 0.1
    
    # pull out the QC 
    qc = data['FparLai_QC']
    # find bit 0
    qc = qc & 1
    
    # generate the masked array
    laim = np.ma.array ( lai, mask=qc )
    
    return laim

import glob

year = 2012
tile = 'h17v03'
files = np.sort(glob.glob('files/data/MCD15A2.A%d*.%s.*.hdf'%(year,tile)))

lai = []
for f in files:
    lai.append(read_modis_lai(f))
    
# force it to be a masked array
lai = np.ma.array(lai)

# plot the data
import pylab as plt

# work out a consistent scaling
lai_max = np.max(lai)

for i,f in enumerate(files):
    plt.figure(figsize=(7,7))
    plt.imshow(lai[i],interpolation='none',vmin=0.,vmax=lai_max*0.75)
    # remember filenames of the form
    # files/data/MCD15A2.A2011185.h09v05.005.2011213154534.hdf'
    file_id = f.split('/')[-1].split('.')[-5][1:]
    print file_id
    # plot a jpg
    plt.title(file_id)
    plt.colorbar()
    plt.savefig('files/images/lai_uk_%s.jpg'%file_id)

# now make a movie ...

import os

cmd = 'convert -delay 100 -loop 0 files/images/lai_uk_*.jpg files/images/lai_uk2.gif'
os.system(cmd)